function [X_mean,X_std,X_ar1] = Calc_quarterly_AR1(Phi_Nz,stationary_prob_Nz,X_Nz)
% computes:
% corr((X(t+1)+X(t+2)+X(t+3))/3,(X(t+4)+X(t+5)+X(t+6))/3)
% E[(X(t+1)+X(t+2)+X(t+3))/3]
% Var[(X(t+1)+X(t+2)+X(t+3))/3]

    % E[X]
    X_mean = sum(stationary_prob_Nz(:).*X_Nz(:));

    % X(k) = E[X(t)*X(t+k)]
    % Compute X(1)+2*X(2)+3*X(3)+2*X(4)+X(5)
    X_cross_terms1 = sum(stationary_prob_Nz(:).*X_Nz(:).*...
        (Phi_Nz*(X_Nz(:) + Phi_Nz*(2*X_Nz(:) + Phi_Nz*(3*X_Nz(:) + Phi_Nz*(2*X_Nz(:) + Phi_Nz*X_Nz(:)))))));

    % Compute 3X(0) + 4X(1) + 2X(2)
    X_cross_term2 = sum(stationary_prob_Nz(:).*X_Nz(:).*(3*X_Nz(:) + Phi_Nz*(4*X_Nz(:) + Phi_Nz*(2*X_Nz(:)))));

    X_std = sqrt(X_cross_term2 - 9*X_mean^2)/3;

    X_ar1 = (X_cross_terms1/9 - X_mean^2)/(X_std^2);

end